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SYMBOLS  USED 


A(k,X;a,,T) 

Fourier  coefficient  of  f(x.X;f,  T)  used  in 
appendix  A 

A(t,u))    a     deterministic     modulating     function, 
which  is  unity  for  a  stationary  process 

D  a  parameter  used  in  Barnett's  parameteri- 

zation of  nonlinear  energy  transfer 
that  represents  the  part  receiving 
energy  from  other  components 

E  total  energy 

/  linear  frequency  =  io/2ir 

/o  mean  frequency 

g  acceleration  due  to  gravity 

G  source  function;    also  used   in  Barnett's 

parameterization  to  represent  the  part 
that  transfers  energy  to  other  com- 
ponents 

h,,  a  counting  symbol  used  in  stationarity 

tests 

i  (-1)"2    in   complex    representation,    also 

used  as  an  index  counter 

/  index  counter 

k  wave    number    vector    with    magnitude 

k   =    |k|;   k  also  is  used  as  an  index 
counter 

K  counter  for  data  points 

index  counter 

L  counter  for  data  points 

M  counter  for  data  points 

N  counter  for  data  points 

Ni  N(k,)  action  density  at  wave  number  k, 

used  in  the  Boltzmann  integral  in  ap- 
pendix C 


p  index  counter 

P  pressure    on    the    free    surface    resulting 

from  air  motion;  also  used  as  a  counter 
for  data  points 

r  position  vector 

R„  nth  order  covariance  function 

S„  nth  order  spectrum;  S  without  subscript 

n  represents  a  second-order  spectrum 

T  total  time;  also  representing  a  coupling 

function  in  the  Boltzmann  integral  in 
appendix  C 

x  space  vector  in  Cartesian  coordinates 

X  a  space  parameter  used  in  appendix  A; 

also    used    as    a    Fourier    transform, 
X(w),  for  data  f r(f) 

y  Cartesian  coordinate 

Y  complex  random  function  for  surface  f 

used  in  the  Fourier-Stieltjes  represen- 
tation 

z  Cartesian     coordinate,     directed     verti- 

cally upward 

f  free  surface  displacement 

8  wave  direction 

7T  3.14159  .... 

q  density  of  water 

t  time 

0  velocity  potential 

$  complex  random  function  for  4>  used  in 

the  Fourier-Stieltjes  representation 

co  radian  wave  frequency 
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SPECTRAL  GROWTH  AND  NONLINEAR  CHARACTERISTICS 
OF  WIND  WAVES  IN  LAKE  ONTARIO 


Paul  Chi  Liu 


ABSTRACT.  Recent  studies  have  shown  that  the  growth  processes  of  wind  waves  are  primarily 
associated  with  the  nonlinear  energy  flux  due  to  wave-wave  interactions.  A  detailed  empirical  examina- 
tion of  these  interactions  uses  calculated  unispectra,  bispectra,  and  trispectra  of  continuously  recorded 
wave  data  during  three  episodes  of  growing  waves.  While  the  unispectrum  provides  information  on  the 
energy  content  of  the  frequency  components,  the  bispectrum  and  trispectrum  generally  provide  informa- 
tion on  the  interactive  relations  between  two-  and  three-frequency  components  respectively.  These 
higher-order  interactive  relations  can  be  considered  characterizations  of  nonlinear  interactions.  The 
results  indicate  that  the  peak-energy  frequency  transfers  more  energy  to  the  lower  frequency  components 
than  to  the  higher  ones,  which  is  confirmation  that  unispectral  peaks  shift  progressively  toward  lower  fre- 
quencies during  wave  growth. 


1.  INTRODUCTION 


The  study  of  wind-generated  waves  has  experi- 
enced significant  and  extensive  development  during 
the  last  20  years  or  so.  Ursell's  (1956)  review,  in 
which  he  found  the  state  of  our  knowledge  of  wind- 
wave  generation  profoundly  unsatisfactory  with  re- 
spect to  both  theory  and  observation,  has  often  been 
credited  with  providing  the  major  stimulation  for 
modern  studies.  After  two  decades  and  vast  theoreti- 
cal and  experimental  efforts,  a  recent  review  by  Bar- 
nett  and  Kenyon  (1975)  observed  that 

At  the  time  Ursell  reviewed  the  field,  the 
body  of  theoretical  work  exceeded  that  of 
the  experimental  work,  but  both  were  in  an 
unsatisfactory  state.  Today  the  same  ratio 
holds  in  that  the  theoretical  ideas  are  still 
ahead  of  the  experimental  testing.  In  par- 
ticular, field  observations  relevant  to  wave 
generation  and  dissipation  in  the  oceans  in 
1955  were  nearly  nonexistent.  Today  they 
are  simply  very  scarce. 

Theory  often  leads  experiment  in  science.  In  the 
study  of  wind  waves,  however,  the  difference  is  so 
large  that  experiment  cannot  interact  with  theory 
very  efficiently.  Theoretical  analyses  result  in  mathe- 
matical complexity,  which  prevents  exact  and  practi- 
cable solutions;  at  the  same  time  a  lack  of  basic 
knowledge  of  the  actual  wave  processes  limits  the 
ability  to  develop  new  models.  This  report  presents  a 
detailed  empirical  examination  of  the  temporal  evo- 
lution processes  of  the  energy  spectrum  of  wind- 
generated  surface  waves.  Our  measurements  and 
analyses,  presented  here,  help  satisfy  the  need  for 
basic  information  about  wave  processes. 

Because  several  recent  studies  (Hasselmann  et 
al.,  1973;  Longuet-Higgins,  1976;  Fox,  1976)  stress 


the  importance  of  nonlinear  energy  transfer  in  the 
wave-growth  processes,  we  have  examined  these  im- 
plications empirically.  Since  the  nonlinearity  comes 
from  the  higher-order  terms  in  the  equations  of 
motion,  the  first  step  in  evaluating  the  nonlinearity 
of  a  set  of  wave  data  recorded  from  a  single  station  is 
to  analyze  the  higher-order  moments  of  the  process. 
In  practice,  this  is  equivalent  to  performing  bispec- 
tral  and  trispectral  analyses  of  the  data. 

A  bispectrum  is  the  two-dimensional  Fourier 
transform  of  the  third-order  covariance  function  of 
the  data.  A  trispectrum  is  the  three-dimensional 
Fourier  transform  of  the  fourth-order  covariance 
function  of  the  data.  As  two-  and  three-dimensional 
Fourier  transforms  of  the  corresponding  covariance 
functions  are  very  cumbersome,  we  have  found  a 
way  to  calculate  bispectra  and  trispectra  by  using  a 
fast  Fourier  transform  algorithm  directly  on  the  wave 
data.  (The  details  will  be  discussed  in  section  5.) 
Physically,  just  as  the  unispectra  provide  informa- 
tion on  the  energy  content  of  the  frequency  compo- 
nents, the  bispectra  and  trispectra  provide  informa- 
tion on  the  interactive  relations  between  two-fre- 
quency components  and  between  three-frequency 
components,  respectively.  We  consider  these  higher- 
order  relations  to  be  estimates  or  characterizations  of 
nonlinear  interactions. 

Hasselmann,  Munk,  and  McDonald  (1962)  and 
Garrett  (1970)  have  demonstrated  that  calculations 
of  observed  bispectra  of  ocean  waves  correlate  rea- 
sonably well  with  theoretically  derived  bispectra. 
These  studies,  however,  are  not  extensive  enough  to 
provide  much  insight  into  the  detailed  behavior  of 
the  physical  processes. 

Using  continuous  wave  data  recorded  in  Lake 
Ontario,  this  report  examines,  identifies,  and  re- 
solves the  temporal  characteristics  of  linear  and  non- 
linear interactions  during  wave  growth.  Our  demon- 
stration that  the  mean  energy  of  the  waves  is  closely 


related  to  the  unispectra,  bispectra,  and  trispectra  of 
the  data  provides  some  analytical  background  to  our 
basically  intuitive  approach. 

2.  THEORETICAL  CONSIDERATIONS 

All  theoretical  studies  of  surface  waves  start 
from  basic  equations  of  fluid  dynamics  with  varied 
idealizations  or  assumptions  for  obtaining  the  solu- 
tions. We  assume  irrotational  motion  of  a  hori- 
zontally unbounded  incompressible  fluid  with  infi- 
nite depth  and  a  free  surface  at  z  =  f  (x,  f),  where  x 
=  (x,  y)  and  z  are  Cartesian  coordinates  and  the  z- 
axis  is  directed  vertically  upward;  then  a  velocity 
potential  0  (x,  z,  f)  exists  and  the  motion  is  governed 
by  the  Laplace  equation 


V20  =  0,  forz  <  f . 


(1) 


Further  neglect  ng  the  surface  tension,  the  kinetic  and 
dynamic  boundary  conditions  at  the  free  surface 
z  =  f  are  given,  respectively,  by 


_^+vr-v0  =  o 


dt 


and 


£+*+>>" 


(2) 


(3) 


where  p  represents  the  pressure  on  the  free  surface  re- 
sulting from  the  air  motion,  q  the  water  density,  and 
g  the  acceleration  due  to  gravity.  Now  if  we  assume 
that  the  random  surface  displacement  is  statistically 
stationary  with  respect  to  both  space  and  time,  we 
can  use  the  Fourier-Stieltjes  representation  (Phillips, 
1966) as 


f(x,f) 


dY(k,u)e 


i(k-x-o)() 


(4) 


Under  the  same  assumptions  and  the  further  assump- 
tion that  <!>—■  0  as  z  —  -oo,  the  solution  to  (1)  can  be 
represented  by 


0(x,z,f)  =  Jki(od$(k/co)ekz-''(k-x-"t). 


(5) 


Here  cM>(k,co)  and  dY(k,u>)  are  complex  random  func- 
tions of  the  horizontal  wave  number  vector  k  and  fre- 
quency co  with  k  =  |k|.  Analogously,  the  pressure 
on  the  surface  resulting  from  air  motion  can  be  repre- 
sented by 


p(x,t)  =  J  J    dP(k/u)e-;(k,x-"rf). 


(6) 


The  representations  (4),  (5),  and  (6)  are  quite  general 
and  appropriate  for  studying  stationary  random 
wave  processes.  Most  of  the  studies  in  the  literature 
consider  processes  that  are  stationary  only  with  re- 
spect to  space  and  thus  reduce  the  equations  of 
motion  into  differential  equations  with  time  deriva- 
tives only.  Examples  are  Phillips  (1960),  Hasselmann 


(1962),  and  Benney  (1962).  Assuming  the  process  is 
stationary  with  respect  to  both  space  and  time,  how- 
ever, further  reduces  the  equations  of  motion  into 
algebraic  operations.  Substituting  (4),  (5),  and  (6) 
into  (2)  and  (3);  expanding  the  factor  ek^  in  (5)  as  a 
power  series;  and  using  (4)  again  in  the  series,  we 
obtain 

-iWY(k,u)  -  /«M>(k,co) 

-  fk,  \u,k' 2d$(k'  ,«'  )dY(k  -  k',« -  co' ) 

-\v\JvA^d^k',<,')dY{k"^") 
x  dY(k-k'  -k",o>-aj'  -co") 

+  Jk,Ju,  -k'  •  (k-k')eM>(k',co') 
x  dY(k-k',co-a>') 

+  \k,\Jy,\ia„-k>.k"k>d*{k',o>') 

xdY(k",w")dY(k-k'  -k",  go -co'  -co") 

=  0, 
and 


-  iW*(k,w)  +  gdY(k,co) 

+  $k,  Jw,{-xV  k'd$(k',w'  )dY(k  -  k\  co  -  co' ) 

+j[k'\k-k'\-k'  .(k-k')]d$(k,o}) 
x  d$(k-k' ,  co-  co')} 

+  k>L'k"L"  [-*«/^y^(k',«/)dY(k",w") 

X  dY(k  -  k'  -  k",  co  -  co'  -  co"  ) 
+  *(*:'  +  k")(k'k"  -k'.k")rf$(k',co') 
X  c/$(k",co")^Y(k-k'  -k",co-co'  -co")"| 
=  --dP(k,co) 


(7) 


(8) 


Equations  (7)  and  (8),  retaining  terms  to  the  third 
order,  can  be  solved  for  dY(k,w)  and  d$(k,co)  if  the 
representation  dP{k,u>)  is  given.  From  the  point  of 
view  of  the  empirical  study,  the  higher  order  terms  in 
(7)  and  (8)  clearly  point  to  the  need  for  study  of 
higher  order  spectra.  Most  of  the  nonlinear  studies  in 
the  literature  assume  the  processes  are  undisturbed 
by  air  motion;  hence  dP(k,co)  =  0  in  (8).  Linearized 
analyses  by  Phillips  (1957)  and  Miles  (1957)  con- 
sidered dP(k,o))  to  be  primarily  associated  with  the 
turbulent  wind  field  and  to  consist  of  two  parts:  the 
part  produced  by  the  turbulent  eddies  in  the  wind 
and  in  random  phase  with  the  wave  field,  and  the 
part  directly  induced  by  and  phase-locked  with  the 
wave  field.  A  combination  of  the  two  mechanisms 
(Phillips,    1966)    shows    that    the   growth    of   wave 


energy  is  initially  linear  with  time  under  the  turbulent 
pressure  alone  and  subsequently  becomes  expo- 
nential because  of  the  induced  pressure  acting  on  the 
growing  waves.  A  more  realistic  representation  for 
dP(k,<j))  is  not  available. 

The  Fourier-Stieltjes  transforms  dY(k,cc)  and 
cM>(k,cd),  solvable  from  (7)  and  (8),  are  closely  related 
to  the  distribution  of  energy  per  unit  among  the  com- 
ponents of  the  wave  field.  The  mean  energy  of  the 
waves  per  unit  projected  surface  area  is 


F  =  F  ■  ,  +  F, 

"-         '-potential     ' 


kinetic 


=  ^egf2+flLe(VcA)2rfz,      (9) 

where  the  bars  denote  ensemble  averages.  It  can  be 
shown  (e.g.,  Phillips,  1961)  that 

hence,  substituting  (4)  and  (5)  into  (10),  we  obtain 
E   =yclkIJk.J<1){g^(k',«')dY(k-k',w-w') 


+  i'w'dY(k',u')d*(k-k',u-w') 
+  LA  „zo>'fc"rfY(k',a>')rf*(k",w") 

J  k     J  a; 

X  dY(k  -  k'  -  k",w  -co'  -  cu"  ) 


...], 


(11) 


The  relation  between  d<l>(k,co)  and  dY(k,co),  correct  to 
the  first  order,  is  given  by 


d<S>(k,co) 


~7^Y(k,co). 
k 


(12) 


Furthermore,  since  we  assume  stationarity  with  re- 
spect to  both  time  and  space,  we  expect  non-zero 
contributions  to  the  integral  in  (11)  to  occur  only 
when  both  k  =  0  and  co  =  0;  thus  (11)  becomes 


=  yeJkju[(2g)dY(k,U)dY*(k,u) 


+  J„  J    uw1rfY(k/«)dY(k1,«i)dY*(k  +  ka,co  +  Wl) 

ki,  J  o)i 


WCOi 


^K)M^gdY^)dY{v^] 


x  dY(k2,co2)dY*(k  +  k2  +  k2,co  +  Ul  +  co3) 

+  ...],  (13) 

where  dY*(k,co)  is  the  complex  conjugate  of  dY(k,co). 
From  appendix  A,  we  see  that  the  averages 
dY(k,co)dY(ki,coi)  .  .  .  dY*(kn,wn)  relate  to  the  nth 
order  spectrum.  Thus,  to  study  the  complete  energy 
distribution,  it  is  necessary  to  examine  the  higher 
order  spectra. 


3.  OBSERVATIONAL  CONSIDERATIONS 

The  results  presented  in  section  2  and  appendix 
A  are  derived  for  the  characterization  of  the  whole 
lake  surface.  To  correlate  the  results  with  actual 
wave  observations,  usually  made  from  a  wave  gage 
installed  at  a  single  selected  location,  the  equations 
are  integrated  and  normalized  over  all  directions  and 
wave  numbers;  this  explicitly  depends  on  time  and 
wave  frequency  only.  We  can  rewrite  equations  A. 8 
and  A. 9  as 


dY(w,)rfY(a>2)  .  .  .  rfY(co„_1)rfY*(a)fI) 


du)idu)2  ■ 


■>n-l 


.  .  \R„(t1m  .  .  .  ,7-„_!)e 


n-\ 

II 
1=1 


(2tt)"-i 
X  dri  .  .  .  dr„_i 

=  Sn(w!,co2,   .   .   .  ,«„_].), 

if  co„  =  coi  +  co2  +  .  .  .   +  wn_i; 

and 

dY{^)dY{u>2)  .  .  .  dYju^dY*^) 
du)ido)2  ■  ■  ■  ^<jJ„_: 

=?  0,  if  co„  it  ah  +  co2  +  .  .  .  +  &)„_!■ 
Equation  (13)  now  becomes 
1    r2gdY(w)dY*(co) 


(14) 


(15) 


+  j    cowj  rfY(co)dY(co1)rfY*(co  +  wx) 


CUCOi 


h  |    \    ^-dY{u)dY{uMY{u2) 


] 


X  rfY*(oo!  +  co2  +  co)  + 

=  yg  \     2gS2(co)dco  +  j    GjaJiS3(co,cJi)^todu)i 


COCOi 


+ 


wi  ^2   2  g 


S4  ( co ,  cOi ,  cu2 )  dudwi  d(ji2 


(16) 


where  S2(co),  S3(co,tOi),  and  S4(co,coi,co2)  are  the  uni- 
spectrum,  bispectrum,  and  trispectrum,  respectively. 
For  a  completely  stationary  process  with  respect 
to  both  time  and  space,  E  is  constant  in  (16).  How- 
ever, our  intent  is  to  examine  the  time  rate  of  change 
of  the  mean  energy  during  a  growing  wave  process 
that  is  in  general  nonstationary.  Because  in  the  non- 
stationary  case  E  changes  with  time,  its  application 
seems  to  be  somewhat  contradictory.  To  pursue  our 
interest  within  the  framework  of  (16),  we  use  a  basic 


assumption  of  local  stationarity  in  which  the  process 
is  assumed  to  be  stationary  and  E  constant  only 
within  a  given  segment  of  data,  a  segment  long 
enough  to  provide  sufficient  degrees  of  freedom  to 
lead  to  reasonable  spectral  estimates.  Although  Eand 
hence  the  spectra  are  constant  within  the  given  seg- 
ment, they  invariably  change  from  segment  to  seg- 
ment. Therefore  if  we  calculate  S2,  S3,  and  S4  over 
consecutive  overlapping  segments  during  wave 
growth,  we  can  indirectly  examine  the  temporal 
growth  behavior  of  the  energy  spectra. 

The  above  scheme  is  demonstrated  graphically 
by  figure  1.  The  discrete  time  series  rj(tn),  tn  =  nAt,  n 
=  1,  2,  .  .  .  ,  N  represents  an  episode  of  growing 
wave  data  NAt  in  total  time  length  with  a  sampling 
rate  of  1/Af.  We  subdivide  r)(tn)  into  M  overlapping 
segments  tm(t(),  t(  =  (At,  I  =  1,  2,  .  .  .,  L,  m  =  1,  2, 
.  .  .,  M.  Each  fm(tp)  is  LAt  in  time  length.  The  data  ( 
=  K  to  t  =  L  of  £m(t()  overlap  the  data  f  =  1  to  i  =  L 
-  K  +  1  of  {•m+i(ff).  Thus  we  have  N  =  (M  -  1)K  + 
L.  While  the  process  rj(tn)  is  not  stationary,  we  as- 
sume each  fm(ff)  is  locally  stationary  and  hence  all 
considerations  we  have  presented  pertaining  to  £"(f) 
can  be  applied  to  fAf(£Af),  1  =  1,2,..  .  .,  L. 

Since  we  represented  f(f)  as  a  continuous 
process, 


where 


To  use  the  discrete  time  series,  we  write 


fA,(f  At)  =  ^/A^'^Y^co), 


(17) 


(18) 


dYAf(co)  =  £  dY(co  +  2-ire/At) 


I 


(19) 


can  be  considered  as  complex  amplitudes  at  frequen- 
cies co  +  2-k(  / At,  in  view  of  the  fact  that  exuAb*  is 
periodic  at  frequency  co  +  2ir/At.  The  relationship 
between  the  continuous  nth  order  spectrum  S„(coi, 
.  .  .,  co„)  and  the  discrete  nth  order  spectrum  S„  Af(coi, 
.  .  .,  as „)  has  been  given  by  Brillinger  and  Rosenblatt 
(1967)  as 

Sn,A((c°l'   •    •    •'  wn) 

=  z,..zsnL+l%± ^  +  ^V(20) 

with 

0J„   =  CO!   +   .    .    .    +  C0„_!, 

and 

«-l 


2*7„  _  V    /        2jr/A 
J"  +    -AT  "  ti  V  '     Af  / 


Thus  using  the  Nyquist  frequency  of  ir/ At  radians 
per  second  or  1/(2  At)  Hz  in  sampling  to  avoid  alias- 
ing problems,  we  can  effectively  represent  the  con- 
tinuous process  with  a  discrete  time  series  in  our  sub- 
sequent analysis. 

Some    further    remarks    concerning    our    sta- 
tionarity assumption  can  be  stated  here: 


^M(t,)-L2 
Cn-1<t|)^-  ■■■£ 
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Figure  1. — Graphical  demonstration  of  data  segments. 


(1)  The  use  of  the  local  stationarity  assumption  is 
not  new;  Kolmogoroff  (1941)  used  local  homogeneity 
in  his  turbulence  studies.  Most  wave  growth  studies 
in  the  literature  use  this  assumption  implicitly.  To 
evaluate  the  approximate  validity  of  this  assumption, 
appendix  B  gives  some  nonparametric  stationarity 
tests  on  the  data  used  in  this  study.  The  results  show 
that  the  assumption  is  generally  within  the  acceptable 
range  of  the  tests  and  can  be  considered  valid  for  our 
purpose. 

(2)  We  are  not  attempting  explicitly  to  calculate  E 
in  (16).  Rather  we  use  (16)  to  show  the  significance  of 
the  spectral  terms  in  connection  with  mean  energy 
and  thus  concentrate  our  interests  on  the  temporal 
growths  of  the  frequency  components  of  S2(co), 
S3(co,a)i),  and  S4(co,coi,u>2)  and  their  physical  implica- 
tions. 

(3)  A  different  approach  to  analyzing  a  non-sta- 
tionary process,  perhaps  more  rigorous  statistically, 
was  developed  by  Priestly  (1965).  He  introduced  a 
representation  of 


f(f)  =  \uA(t,a)dY(a)e- 


(21) 


where  A(t,oi)  is  a  deterministic  modulating  function 
that  approaches  unity  when  the  process  approaches 
stationarity.  This  representation  could  have  been 
used  in  our  analyses.  However,  since  A{t,o>)  is  not  a 
known  function  and  its  application  inevitably  re- 
quires further  assumptions  and  complications,  we 
did  not  use  it.  It  is  of  interest  to  note  that  Priestly's 
model  leads  to  the  conclusion  that  "the  evolutionary 
spectrum  at  each  instant  of  time  may  be  estimated 
from  a  single  realization  of  a  process."  This  is  exactly 
what  the  local  stationarity  assumption  implies. 

(4)  We  define  local  stationarity  here  in  its  literal 
sense.  Silverman  (1957)  introduced  the  concept  of  a 
locally  stationary  random  process  with  a  locally  sta- 
tionary covariance  that  can  be  written  as  the  product 
of  a  stationary  covariance  and  a  nonnegative  func- 
tion. Since  we  allow  our  local  covariance  to  vary 
from  segment  to  segment  and  the  variation  is  gener- 
ally smooth  and  gradual,  it  seems  analogously  pos- 
sible that  a  nonnegative  function  could  depend  on 
time  segments  and  that  a  stationary  covariance  can 
be  deduced  from  our  consecutive  local  covariances. 
As  we  expect  the  property  to  apply  to  third-  and 
fourth-order  covariances  also,  our  assumption  car- 
ries a  different  sense  than  Silverman's  rigorous 
process. 


Figure  2. — Location  of  wind  and  wave  gages  in  Lake  Ontario. 


were  deployed  for  the  IFYGL  programs.  The  wave 
data  in  this  study  were  recorded  from  the  two  Wave- 
riders  designated  as  OS-1  and  OS-2  (figure  2);  the 
corresponding  wind  data  used  here  were  recorded 
from  PDCS  buoy  11. 

The  Waveriders  were  deployed  in  150  m  of 
water  and  freely  moored  to  a  chain  sinker  with  a 
mooring  line  approximately  twice  the  depth  of  the 
water.  The  Waverider  Buoy,  manufactured  by  Data- 
well,  Holland,  is  of  spherical  shape  1  m  in  diameter 
and  weighs  about  100  kg.  It  contains  two  main  com- 
ponents: an  accelerometer  and  a  transmitter.  The  ac- 
celerometer,  mounted  on  a  pendulous  system,  meas- 
ures the  vertical  component  of  acceleration  as  the 
buoy  moves  with  the  waves.  Two  electronic  inte- 
grators in  cascade  then  transform  the  output  into  a 
voltage  that  represents  the  vertical  displacement  of 
the  buoy.  This  voltage  controls  the  frequency  of  an 
audio  oscillator,  which  in  turn  modulates  a  crystal- 
controlled  transmitter  that  transmits  the  signal  by 
telemetry  to  a  shore  receiver.  The  telemetered  wave 
data  were  recorded  continuously  on  analog  magnetic 
tapes.  The  analog  wave  data  tapes  were  subsequently 
processed  through  a  computer  digitization  and  edit- 
ing system  (Liu  and  Robbins,  1974)  to  obtain  final 
data  tapes,  which  were  digitized  at  a  sampling  rate  of 
approximately  three  per  second.  This  sampling  rate  is 
more  than  sufficient  to  avoid  an  aliasing  problem, 
since  the  buoy  response  is  such  as  to  adequately 
damp  waves  having  a  frequency  of  >  1.0  Hz.  On  the 
other  hand  the  Waverider's  frequency  range  is  given 
as  between  0.065  Hz  and  0.50  Hz;  hence,  in  the  ac- 
tual computations  we  use  a  sampling  period  of  Af  = 
2/3  s  to  yield  a  Nyquist  frequency  of  0.75  Hz  for  the 
computed  spectra.  Referring  to  the  scheme  discussed 
in  the  last  section,  we  use  the  following  numbers  in 
the  analysis: 


4.  DATA  ACQUISITION  AND 
PROCESSING 

The  data  used  in  this  study  were  recorded  in 
Lake  Ontario  from  1  April  1972  to  31  March  1973, 
the  International  Field  Year  for  the  Great  Lakes 
(IFYGL).  Seven  Waverider  Buoys  and  a  large  number 
of  Physical  Data  Collection  System  (PDCS)  buoys 


N  = 

15975 

M  = 

64, 

L  = 

1800, 

K  = 

225. 

Consequently  each  selected  episode  is  177.5  minutes 
long.  The  episodes  are  subdivided  into  64  equal  seg- 
ments of  20  minutes  each,  with  an  overlap  of  17.5 
minutes  from  one  segment  to  the  next. 


5.  SPECTRAL  COMPUTATIONS 


In  the  previous  discussion  we  have  been  using 
general  nth  order  spectra.  For  the  actual  applications, 
however,  we  shall  concentrate  only  on  n  =  2,  3,  and 
4  for  the  unispectrum,  bispectrum,  and  trispectrum, 
respectively,  defined  as 


S2(c 


dY(u)dY*(co) 


=  \    R2(T)e-,ulTdi 


(22) 


R4(-T2,T3  ~  T2,Ti  -  T2) 
R*(ti  -  T2,-T2,T3  -  T2) 
R*(ti  ~  T2rT3  -  T2,-72) 
RiW3  -  T2,~T2,Ti  -  T2) 
R*{t3  ~  72,Ti  -  T2,-T2) 
Ri{-Ti,Ti  -  T3,T2  -  T3) 
R^{-T3,T2  -  T3,Ti  -  T3) 
Ro(t2  ~  T3,-T3,Ti  -  Ti) 
Ri(T2  ~  T3,Ti  -  T3,-T3) 
Ra(Ti  ~  Ti,-Ti,T2  -  T3) 
R4(Ti  -  Ti,T2  -  Ti,-T3). 


(30) 


The  corresponding  transforms  of  the  symmetry  rela- 
tions are 


S3(Wi, 


and 


U2) 


dYMdYiwiWiw*  +  oj2) 


S2(o))    =  S2(-oo), 


(31) 


du)idoj2 
=  i    i    Ri(TuT2)e-'{^  +  ^^dTldT2,  (23) 


S4(a>i,co2,u)3)  — 


rfY(a)1)^Y(w2)c/Y(uj3)rfY*(co1  +  Ui2  +  co3) 


do)idu)2du)i 


=     J       J       j      R,(TllT2,Ti)e~'{^T'  +  ^T'  +  UJ'T') 


a)i     OJ2     to 

xdTidr2dT3 


where 


Rz(t)  -  f(f)f(f  +  r), 


(24) 


(25) 


R3(ti,t2)  =  r(t)r(*  +  Tx)f(f  +  t2),      (26) 


and 


R*(TlfTifT3)  -  f(f)f(f  +  Ti)f(r  +  r2)f(f  +  T,).(27) 

The  above  definitions  lead  to  the  following  sym- 
metry relations: 


R2(t)  «K2(-t), 


(28) 


Ri(Ti,T2)    =  R3(t2,ti)  =  R3(-T2,Ti  -  T2) 

=  R3(Ti  -  T2,  -  T2) 

=  R3(-tut2  -  n)  =  R3(t2  -  n,  -  n),     (29) 
and 

RATi,T2lT3)    =  Ri{Ti,T3,T2)   =   Ri{T2,Tx,T3) 
=  Ri{T2lT3,Tl) 

=  R^T3,T-i,T2)   =  Ri{T3,T2,Tl) 
=  Ri{-Tl,T2  -  TUT3  -  7i) 
=  Ri(-TlrT3  -  TUT2  -  Ti) 
=  Ro(t2  -  TU  -  Tx,Ti  -  Ti) 
=  R4(t2  -  Ti,T3  -  Ti,-Ti) 
=  Ra(t3  -  Tlf-TlfT2  -  Ti) 

=  R4(t3  -  ti,t2  -  n,-Ti) 

=  Ri(-T2,T-i  -  T2,Ti  -  T2) 


■>i{U)uO)2) 


=  S3(co2,a;i)  =  S3(co2,-Wi  -co2) 

=  S3(-uh  -  CU2,0>2) 

=  S3(a)i,-C0i  -co2)   =  S3(-OJi 


o)2,03i),      (32) 


and 


S4(coi,w2,co3)   =  S4(aJi,co3,a;2)  =  S4(co2,ah,a>3) 
=  S4(cj2,co3,oji) 

=  S4(co3/u>i,co2)  =  S4(co3,a;2,Wi) 
=  S4(-u>i  -  a>2  -  co3,a)2,a;3) 
=  S4(-wi  -  a>2  -  cj3,w3,w2) 

=  S4(cJ2/-G0i  -  C02  -  C03,O)3) 
=  S4(w2,O)3,-C0i  -  oo2  -  co3) 
=  S4(w3,-a)i  -  a>2  -  oj3,o)2) 

=  S4(c03,C02,-GJi  -  U)2  -  (jj3) 

=  S4(-Wi  -  w2  -  co3,coi,a;3) 
=  S4(-C0i  -  lo2  -  co3,co3,a)i) 
=  S4(co1,-a)1  -cu2  -  C03,U>3) 
=  S4(aJi,O)3,-C0i  -  co2  -  a>3) 
=  S4(aj3,-a)i  -  co2  -  co3,uh) 
=  S4(co3,o)i,-a)i  -  co2  -  co3) 
=  S4(-a)i  -  co2  -  co3,oji,oj2) 

=  S4(-C0i  -  C02  -  C03,C02,O)i) 

=  S4(a)i,-o)i  -  a>2  -  OJ3,C02) 
=  S4(c«Ji,a;2,-coi  -  co2  -  co3) 

=  S4(c02,-Wi  -  C02  -  C03,C0i) 

=  S4(co2,a)i,-cOi  -  a)2  -  oj3). 


(33) 


Because  of  these  symmetries,  we  need  only  to  esti- 
mate the  Sn's,  n  =  2,  3,  and  4,  within  a  fundamental 
region.  The  fundamental  region  for  unispectrum 
S2(co)  is  the  line  segment  0  <  w  <  u)N;  for  bispectrum 
S3(wi,co2)  is  the  triangle  defined  by  0  <  co2  <  ah,  and 
0  <  ooi  <  wN;  and  for  trispectrum  S4(co1,a)2,a;3)  is  the 
tetrahedron  defined  by  0  <  co3  <  w2,  0  <  co2  <  coi, 
and  0  <  o»i  <  ww,  with  ooN  =  2ir/(2At)  representing 
the  Nyquist  frequency. 

Comparing  the  definitions  (22),  (23),  and  (24) 
with  (A. 19)  in  appendix  A,  we  see  that  for  a  finite 
segment  of  a  time  series  the  coefficients  of  its  Fourier 
transform  X(u>)  can  be  used  as  approximations  for  the 
theoretical  values  c/Y(oj)  and  hence  for  the  appro- 
priate products  of  Fourier  coefficients  for  estimating 


corresponding  spectral  densities.  We  used  the  follow- 
ing approach  generally  similar  to  those  discussed  by 
Haubrich  (1965)  and  Hinch  and  Clay  (1968) . 

Starting  with  f(f  At),  2  =  1,  2,  .  .  .,  L,  we  sub- 
divide the  series  into  P  non-overlapping  groups  each 
of  length  K,  such  that  for  p  =  0,  1,  .  .  .,  P  -  1  and  k 
=  1,  2,  .  .  .,  K  we  have 

fp(fcAf)  =  f[(pX  +  /c)Af], 

and  the  K  complex  Fourier  coefficients  for  each  group 
are  given  by 


X>)  =    ^)     EfpO'Ary2^, 


fc  =  1,  2 K.  (34) 

We  can  then  sum  and  obtain  the  average  estimates  by 


wind  conditions  during  these  periods,  recorded  at  10- 
minute  intervals  at  4  m  above  the  lake  surface  on 
PDCS  buoy  11,  are  shown  in  figures  3,  4,  and  5.  The 
group  of  short  and  straight  lines  plotted  on  the  fig- 
ures under  wind  direction  and  wind  speed  indicate 
the  locations  in  time  of  the  64  overlapping  segments 
of  wave  data  analyzed  for  each  episode.  Each  seg- 
ment is  20  minutes  long  and  has  a  17.5-minute  over- 
lap with  the  next  segment.  Although  each  of  the  epi- 
sodes represents  a  growing,  nonstationary  wave 
field,  the  20-minute  segments  are  assumed  to  be 
locally  stationary,  and  thus  we  can  calculate  unispec- 
trum,  bispectrum,  and  trispectrum  for  each  segment 
using  equations  (35),  (36),  and  (37),  respectively. 

The  three  episodes  are  representative  of  growing 
wave  conditions  in  the  Great  Lakes.  Because  we  are 
interested  in  the  early  stage  of  wave  growth,  none  of 
the  episodes  is  under  severe  storm  conditions.  The  9 
August  episode  started  at  a  wind  speed  of  8  m  s"1  and 
increased   to  more   than  11  m  s"1   in  3  hours.   The 


l     L 
S2(oo)  =—     L,  X0(w)Xu(u), 


(35) 


-l 

S3(coi,w2)  =-p      S  Xp(co1)Xp(co2)Xp(coi  +  to2),         (36) 
1      p  =  l 

and 

1       P 
S4(coi,w2,a)3)  —~jS      ^  Xp(co1)Xp(a)2)Xp(aj3) 
1     p=\ 


XXp(coi  +  w2  +  a>3). 


(37) 


Equations  (34)-(37)  represent  the  basic  approach  in 
spectral  computations  used  in  this  report.  The  ap- 
proach is  quite  efficient  and  feasible,  especially  with 
the  fast  Fourier  transform  algorithm  (Cooley  and 
Tukey,  1965)  available  as  a  computing  subroutine.  In 
the  actual  computations,  we  use  P  =  30  for  each 
data  segment  of  L  =  1800  to  obtain  smoothed  spec- 
tral estimates  with  60  degrees  of  freedom  and  a  95 
percent  confidence  interval  between  1.48  and  0.72. 


6.  RESULTS  AND  DISCUSSION 
6.1  The  Episodes 

From  the  ample  supply  of  wave  data  recorded  in 
Lake  Ontario  during  IFYGL,  we  selected  three  epi- 
sodes for  our  study: 

9  August  1972  EST  1215-1515, 

30  September  1972    EST  0540-0840, 
7  October  1972  EST  0000-0300. 

The  August  and  October  episodes  were  recorded 
from  Waverider  OS-2  shown  in  figure  2;  the  Septem- 
ber episode  was  recorded  from  Waverider  OS-1.  The 


9  August  1972 


Figure  3. — Wind  conditions  for  the  episode  of  9  August  1972.  The 
series  of  short  lines  in  the  lower  part  indicates  the  locations  in 
time  of  the  64  overlapping  segments  of  wave  data  analyzed. 


30  September  1972 


7  October  1972 


Figure  4. — Wind  conditions  for  the  episode  of  30  September  1972. 
The  series  of  short  lines  in  the  lower  part  indicates  the  locations 
in  time  of  the  64  overlapping  segments  of  wave  data  analyzed. 


Figure  5. — Wind  conditions  for  the  episode  of  7  October  1972.  The 
series  of  short  lines  in  the  lower  part  indicates  the  locations  in 
time  of  the  64  overlapping  segments  of  wave  data  analyzed. 


westerly  wind  direction  provided  long  and  approxi- 
mately constant  fetches  during  the  episode.  The  30 
September  episode  included  the  passage  of  a  steep 
low-pressure  center  directly  over  eastern  Lake  On- 
tario, and  wind  speeds  increased  from  4  m  s"1  to 
more  than  14  m  s-1  within  an  hour,  while  wind  direc- 
tions were  changing  rapidly  from  south  through  west 
to  north  during  the  early  part  of  the  episode.  The  7 
October  episode  was  milder  with  wind  speed  in- 
creased from  2.5  m  s  "'  to  8.5  m  s  "'  under  south- 
westerly direction.  The  total  wave  energy  under 
these  three  wind  fields  was  growing  steadily  in  each 
of  the  three  cases. 


6.2  The  Unispectra 

We  first  plot  the  computed  unispectral  density 
versus  frequency  versus  time  for  three  episodes  as 
shown  in  figures  6,  7,  and  8.  These  three-dimensional 
perspective  figures  present  a  clear  overview  of  the 


spectral  growth  of  the  episodes.  Several  basic  charac- 
teristics can  be  observed  from  these  figures: 

(1)  The  growth  activity  is  mainly  concentrated 
over  the  low-frequency  side;  the  high-frequency  side 
does  not  change  much  during  the  episode. 

(2)  During  the  growth  period,  the  peaks  of  the 
spectra  tend  to  shift  toward  lower  frequencies. 

(3)  The  growth  rate  varies  for  different  frequency 
components. 

(4)  Once  a  frequency  component  grows  to  be  the 
spectral  peak,  it  reaches  a  relative  equilibrium  and  its 
growth  rate  tends  to  diminish. 

(5)  A  comparison  of  these  local  spectral  growth 
episodes  with  their  respective  wind  conditions  shown 
in  figures  3,  4,  and  5  indicates  that  dominant  growth 
of  the  spectra  happens  during  increasing  wind 
speeds. 

The  above  features  are  generally  known;  similar 
results  were  obtained  in  laboratory  studies,  erg., 
Jacobson  and  Colonell  (1972)  and  more  recently  by 
Wuetal.  (1979). 


Figure  6.  —  Perspective  view  of  the  unispectrum  during  the  episode 
of  9  August  1972. 


Figure  7. — Perspective  view  of  the  unispectrum  during  the  episode 
of  30  September  1972. 


With  respect  to  the  well-known  shifting  of  spec- 
tral peaks  during  wave  growth,  Hasselmann  et  al. 
(1973)  suggested  that  nonlinear  energy  transfer  plays 
a  dominant  role  in  the  process.  They  concluded  that 
the  spectral  form  results  from  the  self-stabilizing 
property  of  the  nonlinear  interactions  that  con- 
tinually readjust  the  energy  distribution  within  the 
spectrum.  We  shall  further  explore  these  im- 
plications. 


6.3  Temporal  Growth  of  Unispectral 
Components 

To  examine  the  time-dependent  behavior  of  the 
frequency  components  individually,  we  first  smooth 
them  in  time  by  hanning  and  then  plot  them  semi- 
logarithmically.  The  components  have  different  time 
dependences,  but  they  can  be  grouped  into  three  fre- 
quency ranges.  The  results  are  shown  in  figure  9  for 
the  9  August  episode,  in  figure  10  for  the  30  Septem- 
ber episode,  and  in  figure  11  for  the  7  October 
episode. 

(1)  Figures  9(a),  10(a),  and  11(a)  show  the  low-fre- 
quency group  that  contains  components  with  fre- 
quencies less  than  all  the  peak-energy  frequencies 
during  the  episode.  The  components  in  this  group 
can  be  characterized  by  their  sensitivity  to  the  wind 
field.  Their  growth  seems  to  follow  directly  with  the 
increases  in  wind  speed.  The  approximate  linearity 
during  the  growth  as  shown  in  the  logarithmic  plot 
indicates  the  exponential  growth  that  Miles'  theory 
predicts. 

(2)  Figures  9(c),  10(c),  and  11(c)  show  the  high-fre- 
quency group  that  contains  components  with  fre- 
quencies beyond  the  peak-energy  frequencies  during 
the  episode.  The  components  in  this  group  represent 
the  portion  of  the  spectrum  that  is  usually  considered 


Figure  8. — Perspective  view  of  the  unispectrum  during  the  episode 
of  7  October  1972. 


to  be  the  equilibrium  range.  They  tend  to  parallel 
each  other  in  time  and  are  insensitive  to  either  in- 
creasing wind  speed  or  time  duration.  Some  anoma- 
lous behavior  shown  in  figure  10(c)  during  the  early 
stage  may  be  due  to  the  fact  that  wind  direction  dur- 
ing this  short  time  interval  was  changing  con- 
tinuously and  thus  further  complicated  the  whole 
picture. 

(3)  Figures  9(b),  10(b),  and  11(b)  show  the  compo- 
nents, lying  between  groups  1  and  2,  that  have  been 
peak  frequency  during  the  episode.  The  behavior  of 
components  in  this  group  is  complicated  since  they 
have  mixed  properties  of  both  of  the  first  two  groups. 
This  group  contains  a  large  part  of  the  total  spectral 
energy. 
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Figure  9a. — Growth  of  low-frequency  unispectral  components  during  the  episode  of  9  August  1972. 
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Figure  9b.  — Growth  of  middle-frequency  unispectral  components  during  the  episode  of  9  August  1972. 
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Figure  9c. — Growth  of  high-frequency  unispectral  components  during  the  episode  of  9  August  1972. 


6.4  The  Source  Function- 
Empirical  and  Theoretical 

One  of  the  main  problems  in  developing  wave 
growth  models  is  determining  the  source  function  G 
in  the  energy  balance  equation  (Gelci  et  al.,  1956; 
Hasselmann,  1960): 


d_E 
dt 


Vv£  +  K 


Vk£  =  G, 


(38) 


where  £  is  the  energy  density,  the  second  and  third 
terms  on  the  left  side  represent  the  advection  and  re- 
fraction of  the  processes,  respectively;  the  dots  de- 
note the  time  derivatives.  The  source  function  G  is 
the  total  rate  of  change  of  the  energy  representing  all 
processes  of  generation,  dissipation,  and  linear  and 
nonlinear  interactions.  Most  wave  studies  in  recent 
years  have  aimed  at  assessing  some  aspects  of  the 
source  function.  Perhaps  the  most  significant  one  is 
the  Joint  North  Sea  Wave  Project  (JONSWAP)  in 
which  Hasselmann  et  al.  (1973)  found  that  the  non- 
linear energy  transfer,  computed  from  integrating  the 
theoretical  Boltzmann  integral  (appendix  C),  is  im- 
portant for  explaining  the  form  of  the  wave  spectrum 
during  wave  growth. 

With  the  available  temporal  functions  of  the  fre- 
quency components,  we  estimated  dE/dt  empirically 
for  each  spectrum  S2(co)  by  fitting  a  smooth  cubic 
spline  curve  through  the  components  and  calculating 


their  time  derivatives.  For  waves  at  a  single  location, 
if  both  the  advection  and  refraction  terms  can  be 
ignored  in  equation  (38),  the  empirical  dS2(w)/dt  also 
represents  the  empirical  source  function.  (Given 
duration  t,  fetch  L,  and  group  velocity  v,  advection 
can  be  neglected  if  t  <  L/v.  This  condition  is  gener- 
ally satisfied  in  this  study  since  duration  is  small  for 
increasing  wind  speed  and  fetch  is  large  during  the 
episodes.) 

Hasselmann  et  al.  (1973)  also  calculated  empiri- 
cal source  functions  by  parameterization  of  (38). 
Their  method  includes  the  advection  term.  However, 
they  give  only  one  average  source  function  for  each 
episode.  Our  method,  on  the  other  hand,  provides 
consecutive  source  functions  throughout  the  episode 
and  thus  more  detailed  information. 

With  these  empirical  source  functions  available, 
the  next  task  is  to  correlate  them  with  the  theoretical 
nonlinear  source  functions.  The  numerical  integra- 
tion of  the  Boltzmann  integral  by  Sell  and  Hassel- 
mann (1972)  requires  30  minutes  of  computer  time 
for  each  spectrum  and,  at  this  stage,  is  rather  imprac- 
tical. In  appendix  C  we  show  that  Barnett's  (1966) 
parameterization  of  Hasselmann's  earlier  calculations 
(1963a, b)  leads  to  essentially  the  same  conclusions  as 
those  obtained  by  Sell  and  Hasselmann.  Therefore 
we  choose  to  use  the  low-frequency  part  of  the  results 
calculated  from  Barnett's  parameterized  scheme  to 
represent  the  theoretical  nonlinear  source  function  in 
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Figure  10a. — Growth  of  low-frequency  unispectral  components  during  the  episode  of  30  September  1972. 
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Figure  10b.— Growth  of  middle-frequency  unispectral  components  during  the  episode  of  30  September 

1972. 
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Figure  10c.  — Growth  of  high-frequency  unispectral  components  during  the  episode  of  30  September  1972. 
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Figure  11a. — Growth  of  low-frequency  unispectral  components  during  the  episode  of  7  October  1972. 
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Figure  lib. — Growth  of  middle-frequency  unispectral  components  during  the  episode  of  7  October  1972. 
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Figure  lie. — Growth  of  high-frequency  unispectral  components  during  the  episode  of  7  October  1972. 
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Figure  12. — Unispectra  ( — )  and  empirical  ( — )  and  Barnett's  (-■■■) 
source  functions  during  the  episode  of  9  August  1972.  The  num- 
bers 2.5,  32.5,  .  .  .,  and  152.5  are  minutes  from  beginning  of 
episode. 
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Figure  13. — Unispectra  ( — )  and  empirical  ( — )  and  Barnett's  (- — -) 
source  functions  during  the  episode  of  30  September  1972.  The 
numbers  2.5,  32.5,  .  .  .,  and  152.5  are  minutes  from  beginning 
of  episode. 


our  subsequent  comparisons.  (More  recently,  Wu  et 
al.  (1979)  also  used  Barnett's  parametric  equation  in 
calculating  nonlinear  energy  transfer  and  found  good 
comparisons  in  the  low  and  intermediate  frequency 
region  of  a  wave  spectrum.)  We  ignore  the  high  fre- 
quency part  of  the  parameterized  results  because  it 
was  based  on  Neumann's  (1953)  assumption  that 
energy  density  is  proportional  to  aT6  at  the  high  fre- 
quency side  rather  than  to  co"5,  which  most  recent 
studies  have  confirmed  as  correct.  The  results  for  the 
three  episodes  are  shown  in  figures  12,  13,  and  14. 
Because  of  gradual  changes  in  the  process,  we  present 
only  six  spectra  and  source  functions  for  each  episode 
at  30  minutes  apart  to  show  the  essential  features.  In 
the  figures  we  plot  the  logarithm  of  unispectra  log 
S2(cj),  and  the  empirical  and  calculated  d[log 
S2(u>)]/dt  all  with  respect  to  the  frequency  of  cc/2-k. 
The  distinction  between  the  empirical  and  calculated 
"theoretical"  source  functions  is  quite  evident.  The 
theoretical    source   function   consistently   has   large 


positive  lobes  at  frequencies  below  the  peak-energy 
frequency.  The  empirical  source  function  develops 
similar  positive  lobes  only  toward  latter  stages  of  the 
episode  when  the  growth  activity  is  intense.  At  the 
beginning  of  the  episode  when  the  growth  is  gener- 
ally slow  or  when  the  waves  are  well  developed  and 
further  growth  is  slowing  down,  the  empirical  source 
function  tends  to  be  much  less  pronounced  in  its 
lobes,  positive  and  negative. 

This  behavior  of  the  empirical  and  theoretical 
source  functions  is  not  particularly  surprising  since 
the  theoretical  results  are  dependent  on  the  shape  of 
the  spectrum  whereas  the  empirical  results  are  de- 
pendent on  temporal  changes  in  the  spectrum.  This 
also  provides  a  possible  explanation  for  the  domi- 
nance of  nonlinear  interactions  during  growth.  If  we 
assume  that  the  empirical  results  represent  the  com- 
plete source  function  and  the  theoretical  results  only 
the  nonlinear  interactions,  their  differences  then  rep- 
resent the  other  processes,  such  as  wave  breaking  and 


15 


2.5" 

4 

- 

| 

- 

N 

X 

II 
II 

•^ 

II 

*X  3 

F 

u 

J  ? 

1      '/ 

1       1 

\ 

</) 

o> 

i      J1 

- 

o    , 

_l 
0 

I    /    | 

1    /      I 
1    /       \ 

- 
■  / 

1          1          1          1 

i     i     i 

X  3- 

E 


32.5" 
1 

H 
J I 

^r^" 

■    lYrU* 

i 

0     01020304050607 
Frequency  w/2n  Hz 


0     01020304050607 
Frequency  u/2nHz 


Figure  14.—  Unispectra  ( — )  and  empirical  ( — )  and  Barnett's  ( ) 

source  functions  during  the  episode  of  7  October  1972.  The 
numbers  2.5,  32.5,  .  .  .,  and  152.5  are  minutes  from  beginning 
of  episode. 


dissipation.  At  the  beginning  of  the  episode,  when 
growth  is  generally  slow  or  when  the  waves  are  well 
developed  and  growth  is  also  slowing  down,  the  non- 
linear interactions  must  be  balanced  by  some  equally 
important  dissipation  process  in  order  to  yield  the 
smaller  total  source  function.  During  rapid  growth, 
on  the  other  hand,  both  the  nonlinear  interactions 
and  the  total  source  functions  have  large  positive 
peaks,  which  implies  that  dissipation  is  less  signifi- 
cant and  that  most  of  the  growth  is  through  non- 
linear interaction. 


6.5  The  Bispectra 

There  are  three  published  bispectral  studies  on 
wind-generated  waves:  Hasselmann,  Munk,  and 
McDonald  (1962),  and  Garrett  (1970)  compared  four 
bispectra  computed  from  actual  wave  data  with 
those   theoretically  derived  and  found  satisfactory 


agreement,  and  Houmb  (1974)  presented  three 
examples  of  computed  wave  bispectra.  In  the  present 
study  we  are  interested  in  the  higher  order  interac- 
tions during  wave  growth.  As  a  first  step  we  have 
computed  64  consecutive  bispectra  for  each  of  the 
three  growth  episodes.  Because  the  temporal  varia- 
tions are  gradual,  we  present  only  six  bispectra  from 
each  episode  at  30  minutes  apart  as  representative  of 
the  whole  process.  Figures  15,  16,  and  17  show  the 
sample  bispectra  for  episodes  9  August,  30  Septem- 
ber, and  7  October. 

Since  the  bispectrum  S3(cdi,to2),  computed  from 
equation  (36),  is  complex,  it  is  convenient  to  express 
it  as  a  bispectral  amplitude  rather  than  as  co-bi- 
spectra  and  quadrature-bispectra  separately.  Figures 
15,  16,  and  17  show  the  contours  of  the  logarithms  of 
bispectral  amplitude  plotted  in  the  fundamental  re- 
gion of  0  <  u>2  ^  w3  and  0  <  «i  <  cow. 

Several  main  features  can  be  seen  in  these 
figures: 

(1)  The  bispectrum  has  a  hill  at  (cop,gop),  where  cjp  is 
the  frequency  of  the  spectral  peak  in  the  correspond- 
ing unispectrum. 

(2)  Each  bispectrum  also  shows  two  ridges  ap- 
proximately parallel  to  the  two  frequency  axes  and 
sloping  down  from  the  bispectral  hill  toward  both 
higher  and  lower  frequencies. 

(3)  The  magnitude  of  the  bispectral  amplitude  in- 
creases and  the  hills  and  ridges  migrate  toward  lower 
frequencies  during  wave  growth. 

(4)  These  observed  characteristics  are  qualitatively 
similar  to  those  depicted  by  Hasselmann,  Munk,  and 
McDonald's  (1962)  deep-water  wave  co-bispectrum. 

It  is  not  readily  discernible  that  these  features  re- 
sult from  a  nonlinear  process.  However,  from  the 
qualitative  resemblance  to  Hasselmann,  Munk,  and 
McDonald's  results  and  the  fact  that  their  results 
agree  with  those  derived  from  a  theoretical  nonlinear 
process,  we  may  conjecture  that  these  features  are 
characteristic  of  a  weakly  nonlinear  wave  growth 
process  and  hence  proceed  to  examine  their  temporal 
behavior. 


6.6  Temporal  Growth  of 
Bispectral  Components 

Since  the  predominant  bispectral  interactions 
occur  when  the  frequency  of  the  spectral  peak  in  the 
corresponding  unispectrum  interacts  with  itself  and 
these  interactions  migrate  toward  lower  frequencies 
during  wave  growth,  it  is  sufficient  that  we  examine 
the  time-dependent  behavior  of  those  bispectral  com- 
ponents on  the  45°  line  whose  frequencies  interact 
with  themselves.  We  proceed  in  a  manner  similar  to 
that  for  unispectral  components  by  temporally 
smoothing  them  with  hanning  and  then  plotting  them 
semi-logarithmically  as  a  function  of  time.  The  re- 
sults are  quite  similar  to  those  unispectral  compo- 
nents. We  again  separate  them  into  three  groups  for 
the  three  episodes  (figures  18,  19,  and  20). 


16 


0.3 


0.2 


0.1 


0 


i    i    i    i    i    !    i    i    i    i    i    i    i    i — i — r~ r 


62.5 


'fit-  •!•    fit >     '      ■<■»■- 


0  0.1  0.2  0.3  0.4  0.5         0.6         0.7 

Frequency  w/2ir  (Hz) 


0.3 

i    i    i    i    i    i    :    i   i 

1     !     1     1     1     1     i 

I     1     1     i     1     I     1     1     1     1     1     !     ! 

g                 92.5: 

0? 

0.1 

Jfi^  ^  '^r        ■ 

n 

0.3 


0.2 


0.1 


I     I     !     I     I     I     I     I 


,'.,■,.).,    .S^;^. 


03 

i   i   i   i   i   i   i   i   i 

III; 

i     l     l     l     I     I     I     I 

152.5: 

;■■             ■■*■■    t 

0.2 

0.1 

.  .' 

,_  —        v 

' 

n 

>'-':'T      f     '■:    '  i'      f     ■  ;      ■*'- 

0  0.1  0.2         0  3         0.4  0.5         0.6         0  7 

Frequency  w/2ir  (Hz) 


Figure  15. — Bispectra  during  the  episode  of  9  August  1972.  The  numbers  2.5,  32.5,  .  .  .,  and  152.5  are 
minutes  from  beginning  of  episode. 


Several  points  should  be  noted  here: 

(1)  The  figures  are  plotted  with  log[S3(co,co)]  versus 
linear  time  as  before  but  with  a  smaller  range.  Conse- 
quently the  exponential  growth  rates  shown  in  fig- 
ures 18,  19,  and  20  that  appear  comparable  with 
those  shown  in  figures  9,  10,  and  11  are  actually 
about  four  times  smaller. 

(2)  We  have  seen  that  in  unispectral  growth  the 
spectral  peak  moves  to  a  lower  frequency  during 
wave  growth.  The  bispectra  provide  some  indication 
of  phenomena  that  can  be  attributed  to  this  transfer, 
namely,  that  the  interactions  of  the  peak  frequency 
with  the  next  lower  frequency  grow  consistently 
stronger  relative  to  the  interactions  with  the  next 
higher  frequency  during  the  latter  part  of  the  growth. 
Examples  are  shown  in  figures  21,  22,  and  23.  As 
growth  continues,  the  component  of  the  next  lower 
frequency  eventually  becomes  the  spectral  peak. 

(3)  If  the  bispectral  amplitudes  represent  the  sec- 
ond-order nonlinear  interactions  in  the  wave  growth 
process,  as  we  have  conjectured,  then  our  results 
show  that  the  nonlinear  interactions  should  also  have 
functional  relationships  with  respect  to  time.  This  is 
not    unexpected.    Since    the    theoretical    nonlinear 


source  term  is  a  function  of  the  unispectrum,  it  must 
vary  with  time  if  the  unispectrum  does. 

6.7  The  Trispectra 

Brillinger  and  Rosenblatt  (1967)  presented  one 
set  of  calculated  trispectra  for  daily  sunspot  num- 
bers. Their  work  represents  the  only  trispectral  calcu- 
lations available  in  the  literature.  They  did  not,  how- 
ever, provide  an  interpretation  of  their  results.  Tri- 
spectral analysis  has  remained  relatively  unexplored. 

By  extending  the  calculation  scheme  Haubrich 
(1965)  used  in  bispectral  calculations  one  step  further 
and  using  equation  (37),  we  have  been  able  to  calcu- 
late a  trispectrum  for  each  of  the  64  data  segments  for 
each  episode.  It  is  not  immediately  clear  how  to 
visualize  and  present  the  trispectra  effectively.  It  is 
even  more  cumbersome  to  try  to  present  them  in 
time.  Since  the  fundamental  region  for  the  trispec- 
trum is  the  tetrahedron  defined  by  0  <  oo3  ^  w2,  0  < 
<j02  <  ah,  and  0  <  coi  <  coN,  we  have  the  trispectral 
density  S4(uh,t02,aJ3)  in  addition  to  the  three  fre- 
quency axes  to  form  a  four-dimensional  manifold  im- 
bedded in  three-dimensional  Euclidean  geometry.  By 
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Figure  16. — Bispectra  during  the  episode  of  30  September  1972.  The  numbers  2.5,  32.5,  .  .  .,  and  152.5 
are  minutes  from  beginning  of  episode. 


cutting  the  tetrahedron  along  one  frequency  axis  into 
discrete  layers  parallel  to  the  plane  formed  by  the 
other  two  frequency  axes  and  contouring  the  loga- 
rithm of  the  trispectral  amplitudes  on  these  layers, 
we  can  represent  each  trispectrum  with  a  set  of  tri- 
spectral sheets.  Figures  24,  D.l,  and  D.2  show  the 
calculated  trispectra  for  each  episode;  each  figure 
contains  six  sets  of  trispectra  30  minutes  apart,  corre- 
sponding in  time  to  the  bispectra  shown  in  figures  15, 
16,  and  17. 

Some  features  are  observed  to  be  analogous  to 
the  bispectral  features: 

(1)  The  highest  trispectral  interactions  generally 
occurred  at  (a^,cop,  u)p),  that  is,  where  the  spectral 
peak  frequency  of  the  corresponding  unispectrum  in- 
teracts with  itself  to  the  third  order. 

(2)  The  interactions  taper  off  from  (w„,ajp,wA7) 
toward  higher  and  lower  frequencies  along  the  three 
frequency  directions. 

(3)  The  intensity  of  the  interactions  increases  and 
migrates  toward  lower  frequencies  during  wave 
growth. 

These  features  represent  natural  generalization  of  bi- 


spectral characteristics.  Without  appropriate  theo- 
retical nonlinear  analysis  as  guidance,  it  is  not  prac- 
tical to  identify  these  features  as  a  nonlinear  process. 
However,  because  the  features  appear  to  be  coherent 
with  those  of  bispectra  and  because  our  bispectral 
features  resemble  Hasselmann,  Munk,  and  Mc- 
Donald's (1962)  theoretical  results,  we  may  again 
conjecture  that  the  trispectral  features  we  have  ob- 
served are  characteristic  of  a  weakly  nonlinear  wave 
growth  process.  If  a  third-order  nonlinear  theoretical 
analysis  is  performed,  perhaps  extending  that  of 
Hasselmann,  Munk,  and  McDonald,  we  may  expect 
that  theoretical  trispectra  would  produce  features 
similar  to  those  shown  in  figures  24,  D.l,  and  D.2. 

6.8  Temporal  Growth  of 
Trispectral  Components 

Using  the  same  procedure  we  applied  to  examine 
the  time  dependence  of  the  unispectral  and  bispectral 
components,  we  first  plot  the  logarithm  of  the  tri- 
spectral amplitudes,  whose  peak  energy  frequency 
interacts  with  itself  to  the  third  order  with  respect  to 
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Figure  17. — Bispectra  during  the  episode  of  7  October  1972.  The  numbers  2.5,  32.5,  .  .  .,  and  152.5  are 
minutes  from  beginning  of  episode. 


time.  The  results,  shown  in  figures  25,  26,  and  27,  are 
again  quite  similar  to  the  unispectral  and  bispectral 
components.  Since  we  have  only  the  lower  frequency 
components,  grouping  is  not  necessary.  The  results 
show  again  that  the  third-order  interactions  have 
functional  relationships  with  respect  to  time  during 
the  wave  growth  process. 

In  figures  28,  29,  and  30,  we  concentrate  on  the 
peak-energy  frequency  achieved  at  the  latter  stage  of 
the  episode  by  plotting  the  temporal  growth  of  the 
frequency  interacting  with  itself  to  the  third  order, 
the  frequency  interacting  with  itself  and  the  next 
lower  frequency,  and  the  frequency  interacting  with 
the  second-order  interaction  of  the  next  higher  fre- 
quency. These  figures  show  quite  clearly  that  the 
shifting  of  the  peak-energy  frequency  toward  the 
lower  frequency  is  due  to  the  stronger  interactions  of 
the  peak-energy  frequency  with  the  lower  frequency 
and  that  the  interaction  provides  energy  transfer 
from  the  peak-energy  frequency  to  the  lower  fre- 
quency. This  confirms  that  the  higher-order  non- 
linear interactions  are  responsible  for  the  shifting  of 
spectral  peaks  during  wave  growth. 


6.9  Further  Remarks  on 

Bispectra  and  Trispectra 

What  have  we  gained  from  bispectral  and  tri- 
spectral  analysis?  First,  we  have  developed  methods 
for  computing  and  presenting  bispectra  and  tri- 
spectra from  data.  This  experience  will  be  useful  in 
the  future.  Second,  and  more  important,  we  have 
presented  evidence  that  during  the  later  intensive 
part  of  wave  growth,  the  interactions  of  the  peak- 
energy  frequency  component  with  the  next  lower  fre- 
quency component  grow  consistently  stronger  com- 
pared with  the  interactions  with  the  next  higher  fre- 
quency component,  and  the  next  lower  frequency 
subsequently  becomes  the  peak-energy  frequency. 
This  result,  which  provides  a  clear  explanation  for 
the  well-known  fact  of  shifting  spectral  peaks  toward 
lower  frequencies,  can  be  demonstrated  only  by  bi- 
spectral and  trispectral  analysis. 

Benjamin  and  Feir  (1967)  showed  theoretically 
and  experimentally  that  for  deep-water  waves  with  a 
fundamental  frequency  a>,  because  of  coupling 
through  the  nonlinear  boundary  conditions  at  the 
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Figure  18a.  — Growth  of  low-frequency  bispectral  components  during  the  episode  of  9  August  1972. 


18 


N 


g    14 

u 


3   10 


*  *  * 

+ 
+  + 


**>  +  +  ; 

*  <» 


i  e  +  +  +    + 


*******     ++++++++ ++ + + + ++ + + + + + 

*  +  +  +  ++  +  +  +  ++**  *  *  *  x  X  X  X  X  X  X  Z.  i  X  X  ** 

+  **XXXXXX  **  *******00 

XXXXXXXXX  0 

X    X    X    X  0  0 


+  +  +  +  +  - 


0  I 

0    000000000000 


0  0  0  0  0  0  0 


00000000 


CO 

&> 
o 


<D1/27r|<u2/27r 

Hz 

0.175 

0 

0.200 

x 

0.225 

+ 

0.250 

* 

fil  1 1 1 1 1 1 1 I 1 1 1 1  1 1 1  1 1 1 1 1 1 1 1 1 I Ill 


0 


30  60  90  120 

Minutes  from  beginning  of  episode 


150 


Figure  18b. — Growth  of  middle-frequency  bispectral  components  during  the  episode  of  9  August  1972. 
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Figure  18c. — Growth  of  high-frequency  bispectral  components  during  the  episode  of  9  August  1972. 
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Figure  19a. — Growth  of  low-frequency  bispectral  components  during  the  episode  of  30  September  }972. 
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Figure  19b. — Growth  of  middle-frequency  bispectral  components  during  the  episode  of  30  September 

1972. 
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Figure  19c. — Growth  of  high-frequency  bispectral  components  during  the  episode  of  30  September  1972. 
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Figure  20a. — Growth  of  low-frequency  bispectral  components  during  the  episode  of  7  October  1972. 
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Figure  20b. — Growth  of  middle-frequency  bispectral  components  during  the  episode  of  7  October  1972. 
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Figure  20c. — Growth  of  high-frequency  bispectral  components  during  the  episode  of  7  October  1972. 
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Figure  21. — Growth  of  bispectral  components  during  the  episode  of  9  August  1972.  Above:  0.175  Hz  vs. 
0.150,  0.175,  and  0.200  Hz;  Below:  0.200  Hz  vs.  0.175,  0.200,  and  0.225  Hz. 


free  surface,  energy  is  transferred  to  the  side-band 
frequencies  to  (1  ±  f).  Our  results  appear  to  re- 
semble, to  some  extent,  this  side-band  energy  trans- 
fer theory.  We  may  conjecture  that,  analogously,  the 
stronger  interaction  with  the  lower  side-band  fre- 
quency component  than  with  the  higher  side-band 
frequency  component  is  also  due  to  the  complicated 
couplings  of  the  nonlinear  boundary  conditions  at 
the  free  surface.  Furthermore,  Longuet-Higgins's 
(1976)  recent  theory,  indicating  that  the  transfer  of 
energy  tends  to  reduce  any  symmetry  in  the  spec- 
trum, can  also  be  qualitatively  attributed  to  these  un- 
equal side-band  interactions. 


7.  SUMMARY  AND 

CONCLUDING  REMARKS 

In  this  report  we  set  out  to  empirically  examine 
the  temporal  growth  processes  of  wind-generated 
waves  using  data  recorded  in  Lake  Ontario.  We  have 
examined  unispectra,  bispectra,  and  trispectra  com- 
puted consecutively  for  the  three  selected  episodes. 
The  scheme  used  for  our  study  is  based  on  the  as- 
sumed property  of  local  stationarity.  This  assump- 
tion allows  us  to  apply  successfully  the  analysis 
method  developed  for  stationary  processes  to  study 
nonstationary  wave  growth  processes. 
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Figure  22. — Growth  of  bispectral  components  during  the  episode  of  30  September  1972.  Above:  0.175  Hz 
vs.  0.150,  0.175,  and  0.200  Hz;  Below:  0.200  Hz  vs.  0.175,  0.200,  and  0.225  Hz. 


Results  can  be  summarized  as  follows: 

(1)  Wind  waves  grow  fastest  during  increasing 
wind  speeds. 

(2)  Each  unispectral  component  has  a  different 
temporal  growth  rate.  The  components  can  be 
grouped  into  three  frequency  groups:  the  low-fre- 
quency components,  which  consistently  grow  under 
increasing  wind  speeds;  the  high-frequency  compo- 
nents, which  are  insensitive  to  wind  speed  and  time; 
and  the  mid-frequency  components,  which  contain 
those  attaining  spectral  peaks  during  the  episode  and 
which  contain  mixed  properties  of  both  high-  and 
low-frequency  component  groups. 


(3)  Bispectral  and  trispectral  components  generally 
demonstrate  temporal  growth  behavior  similar  to 
that  of  unispectral  components  with  smaller  growth 
rates. 

(4)  From  comparing  empirical  source  functions 
with  parameterized  theoretical  nonlinear  source 
functions,  we  found  that  at  the  beginning  and  ending 
stages  of  wave  growth  the  growth  rate  is  low  and  the 
empirical  source  function  is  smaller;  the  larger  non- 
linear source  function  under  the  spectral  peak  must 
be  balanced  by  significant  dissipation  processes.  Dur- 
ing intensive  wave  growth,  where  nonlinear  and  em- 
pirical source  functions  are  both  large,  dissipation 
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Figure  23. — Growth  of  bispectral  components  during  the  episode  of  7  October  1972.  Above:  0.225  Hz  vs. 
0.200,  0.225,  and  0.250  Hz;  Below:  0.250  Hz  vs.  0.225,  0.250,  and  0.275  Hz. 


becomes  less  significant. 

(5)  During  the  latter  part  of  wave  growth,  the  in- 
teractions of  the  peak-energy  frequency  component 
with  the  next  lower  frequency  component  grow  con- 
sistently stronger  than  the  interactions  with  the  next 
higher  frequency  component,  and  the  next  lower  fre- 
quency subsequently  becomes  the  peak-energy  fre- 
quency. This  result,  as  demonstrated  clearly  from  an 
examination  of  the  temporal  growth  of  bispectral  and 
trispectral  components,  provides  an  explanation  for 
the  well-known  fact  of  unispectral  peaks  shifting 
toward  lower  frequencies  during  wave  growth. 

These  results  are  consistent  with  our  primary  in- 


terest of  exploring  the  empirical  aspects.  Since  the  de- 
tailed process  of  wave  growth  is  still  far  from  being 
completely  understood,  we  hope  our  results  will  pro- 
vide some  insight  that  can  be  useful  for  further 
understanding  of  it. 
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Figure  24a. — Trispectra  during  the  episode  of  9  August  1972,  2.5  minutes  from  beginning  of  episode. 
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Figure  24b.— Trispectra  during  the  episode  of  9  August  1972,  32.5  minutes  from  beginning  of  episode. 
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Figure  24c. — Trispectra  during  the  episode  of  9  August  1972,  62.5  minutes  from  beginning  of  episode. 
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Figure  24d.— Trispectra  during  the  episode  of  9  August  1972,  92.5  minutes  from  beginning  of  episode. 
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Figure  24f .  — Trispectra  during  the  episode  of  9  August  1972,  152.5  minutes  from  beginning  of  episode. 
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Figure  25. — Growth  of  trispectral  components  during  the  episode  of  9  August  1972. 
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Figure  26.  — Growth  of  trispectral  components  during  the  episode  of  30  September  1972. 
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Figure  28b. — Growth  of  trispectral  components  during  the  episode  of  9  August  1972:  Third-order  inter- 
actions among  components  0.200,  0.225,  and  0.250  Hz. 
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Figure  29a. — Growth  of  trispectral  components  during  the  episode  of  30  September  1972:  Third-order 
interactions  among  components  0.175,  0.200,  and  0.225  Hz. 
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Figure  29b. — Growth  of  trispectral  components  during  the  episode  of  30  September  1972:  Third-order 
interactions  among  components  0.200,  0.225,  and  0.250  Hz. 
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Figure  30a.— Growth  of  trispectral  components  during  the  episode  of  7  October  1972:  Third-order  inter- 
actions among  components  0.175,  0.200,  and  0.225  Hz. 
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Figure  30b. — Growth  of  trispectral  components  during  the  episode  of  7  October  1972:  Third-order  inter- 
actions among  components  0.200,  0.225,  and  0.250  Hz. 
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Appendix  A:  Relationships  Between  Higher  Order  Covariances  and  Higher  Order  Spectra 


In  order  to  define  a  general  higher  order  spec- 
trum in  terms  of  the  Fourier  components  of  the  sur- 
face displacement,  we  follow  Batchelor's  (1953)  basic 
approach  by  considering  a  stationary  random  func- 
tion of  both  time  and  space  defined  by 


f(x,f)  for    -X  <  xt :  <  X 
and  -T  <  t  <  T 


r(x,X;f,T)=- 


(A.l) 


0 


for    \xj\  >  X 
and  If I  >  T 


which  satisfies  the  condition  for  a  Fourier  integral  to 
exist.  The  corresponding  Fourier  coefficient  is 

A(k,X;co,T)  =    tt^tJ   \t(x,X;t,T)e+i^™»dxdt. 

(A.2) 

From  Weiner's  generalized  harmonic  analysis,  it  can 
be  shown  that  as  the  parameters  X  and  T  extend  to 
oo,  the  limit  exists  as 

[y(k,c)f;:: 


=  lim  £,$«.  J*.  r'AiKX^Vdk.dkzdhdu 


X-oo 

T-oo 


(A.3) 
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We  can  write 


dY(k,cc)  =  [Y(k,o))] 
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^        y  +  idk:X:      i\     /  „-idut 

which  is  the  inverse  relation  of 


f(x,f)  =  j .  (  dY(k,u)e-'^*-« f>. 


(A.4) 


(A.5) 


Now  we  can  obtain  the  ensemble  average  of  an  nth 
order  product  of  c/Y(k,co)  as 


dYiKaJdYik^)  .  .  .  dY(k„.1,w„_1)dY*(k„,w„ 
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Now  for  stationarity  with  respect  to  both  space  and 
time,  the  right  side  can  be  written  as 
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Since  the  stationary  process  is  independent  of  both  x 
and  r,  the  above  relation  is  identically  zero  unless 


and 


k„  -  ki  -  k2-  ...  -  k„_j  —  0 
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The  last  expression  in  (A.8)  gives  the  general  defini- 
tion of  the  nth  order  energy  spectrum,  where 
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For  observations  at  single  points,  we  may  integrate 
and  normalize  the  above  equations  over  all  directions 
and  wave  numbers  to  obtain 


For  n   =  2,  we  have  the  following  familiar  pair  of 
classic  spectral  analyses: 


and 


S2(co)  =  f-  \R2{r)e-™dT, 

1  2tt       l 


R2(t)=  \S2(o))e™dw. 


(A. 13) 
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Detailed  discussions  corresponding  to  the  formal 
definitions  (A.ll)  and  (A.12),  as  well  as  conditions 
for  existence  and  convergence  of  estimates,  are 
beyond  the  scope  of  our  present  study  and  can  be 
found  in  Brillinger  and  Rosenblatt  (1967). 
In  practical  applications,  we  may  define 


fT(r)  =  f(f),  0  <  t  <  7; 
=  0,  t>  T. 


(A. 15) 


Let  the  Fourier  transform  of  f _(f)  be  X(co);  then  we 
have 


X(w)  =  \^T(t)e-^dt, 


(A.16) 


and 


f7(f)  =  $  uX(a>)eiutdu.  (A. 17) 

Now  for  the  nth  order  covariance  we  have 
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Comparing  (A. 18)  with  integrated  (A. 12),  we  have 
formally 

S„(o>i w„-i) 

=  lim  y  [X(u>a)  .  .  .  Xiu^X*^  +  .  .  .  +  a)„_a)]. 
7^°°  (A.19) 

Thus,  in  practical  applications  we  can  use  the 
product  of  Fourier  transforms  of  fr(f),  with  the  aid  of 
the  fast  Fourier  transform  algorithm,  to  facilitate  the 
calculations  of  higher  order  spectra. 
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Appendix  B:  Testing  for  Stationarity 


In  this  report  we  have  assumed  local  stationarity 
to  study  the  temporal  growth  of  spectral  components 
that  are  in  general  nonstationary.  To  test  the  assump- 
tion of  local  stationarity,  we  adopted  a  testing 
method  similar  to  that  used  by  Bjerkaas  (1976).  The 
scheme  divides  each  data  set  into  P  equal-length  sub- 
segments;  calculates  the  f ,  f2,  f3,  and  f4  for  each 
sub-segment;  and  then  applies  two  nonparametric 
tests,  the  run  test,  and  the  reverse  arrangement  test 
(Bendat  and  Piersol,  1966)  to  the  data  set  to  examine 
its  stationarity. 

In  the  run  test,  the  test  parameter  for  each  inter- 
val is  compared  to  the  median  value  of  the  test 
parameter  for  all  intervals.  If  the  test  parameter  is 
greater  than  or  equal  to  the  median  value,  +1  is  as- 
signed to  the  interval;  otherwise,  -1  is  assigned.  A  se- 
quence of  consecutive  +l's  or  -l's  is  called  a  run, 
and  the  number  of  runs  Nr  in  the  P  intervals  is  deter- 
mined. Nr  gives  an  indication  as  to  whether  or  not  re- 
sults are  independent  random  observations  of  the 
same  random  variable.  The  reverse  arrangement  test 
is  most  useful  in  detecting  monotonic  trends  in  the 
time  series.  Consider  the  test  parameter  A\,  with  i  = 
1,  2,  .  .  .,  P.  If  Ai  >  Aj  for  i  <  j,  the  pair  of  param- 
eters is  called  a  reverse  arrangement. 

The  total  number  of  reverse  arrangements,  Na,  is 
defined  as  follows: 


Let 


Then 


and 


Na 


N   = 


p 

/  =  /  +  ! 


P-l 

L 


Na 


Based  on  P  and  the  level  of  significance,  the  ac- 
ceptance ranges  for  Nr  and  Na  can  be  calculated 
(Bendat  and  Piersol,  1966)  and  used  as  a  basis  to  ac- 
cept or  reject  the  stationarity  assumption. 

Accordingly,  in  our  analysis  eight  separate  tests 
were  performed  for  each  data  set.  The  results,  pre- 
sented in  terms  of  the  percentage  of  data  failing  the 
tests  for  the  three  episodes  studied,  are  shown  in 
table  B.l.  At  a  level  of  significance  a  =  0.10,  the 
overall  failure  rate  was  29  percent;  the  rate  reduced 
to  19  percent  for  a  =  0.02.  Every  data  set  passed  at 
least  two  of  the  eight  tests  performed.  A  total  of  27 
percent  passed  all  eight  tests.  The  parameter  f3  posed 
the  best  results;  £*4,  on  the  other  hand,  posed  the 
worst  results.  Although,  ideally,  absolute  sta- 
tionarity requires  that  the  data  pass  all  the  tests,  we 
find  that  for  a  basically  nonstationary  process  these 
results  generally  tend  to  enhance  the  acceptability  of 
our  assumption  of  local  stationarity.  Thus,  without 
pursuing  more  complicated  physical  and  statistical 
nonstationary  analysis,  we  can  employ  available  sta- 
tionary analysis  methods  to  study  an  otherwise  non- 
stationary  process. 


Table  B.l.  Stationarity  test  results 


Percentage  of  Data  Failed 


Au 
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All  Data 

Level  of 

a  =  0.10 

a 

=  0.02 

a  =  0.10 
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a  =  0.10 

a  =  0.02 
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=  0.02 

Significance 
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45 
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36 
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16 
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0 

6 

0 

3 

0 

4 
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r4 

98 

92 

73 

63 

52 

41 

74 

65 

T 

0 

0 
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0 

2 

0 
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0 
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28 

17 

53 

30 
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33 

44 

27 

f3 
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34 

20 

53 
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Appendix  C:  Application  of  Barnett's  Parameterization  of  a  Nonlinear  Source  Function 

Theoretical  results  have  shown  that  the  non- 
linear energy  transfer  of  gravity  waves  can  be  ex- 
pressed by  the  Boltzmann  integral  as 

|^-4  =  j  .  .  .  jr(k1,k2/k3,k4)[(N3  +  NJNyNz 

-  (N:  +  N2)N3Nii]  5(a)!  +  a>2  -  w3  -  co4) 

x  5(ki  +  k2  -  k3  -  k4)dkjdk2dk3  (C.l) 

where  N,  =  N(k;)  =  S(k,)/co,  with  co,  =  (g/c,)1/2  and  S- 
(k,)  is  the  two-dimensional  wave  spectrum  with  re- 
spect to  wave  numbers  |k,|  =  k\.  The  delta  functions 
express  conditions  for  resonance  between  waves.  The 
coupling  coefficient  T(k1/k2,k3,k4)  is  a  very  compli- 
cated function  and  its  precise  form  does  not  provide 
any  physical  interpretation.  While  the  Boltzmann 
integral  has  been  applied  to  other  fields  of  study, 
Hasselmann  (1962,  1963a, b)  first  derived  equation 
(C.l)  for  gravity  waves.  Other  derivations  have  been 
given  by  Benney  and  Saffman  (1966)  and  Watson 
and  West  (1975).  Recently  Longuet-Higgins  (1976) 
presented  a  simplified  model  and  showed  that 
T(k1,k2,k3,k4)  is  equal  to  47r  when  the  four  wave 
numbers  k1#  k2,  k3,  and  k4  are  nearly  equal.  This  re- 
sult implies  that  the  exchange  of  energy  within  the 
peak  of  the  spectrum  is  of  dominant  importance.  No 
practical  application  of  this  simplified  model  has 
been  developed,  however. 

Attempts  to  evaluate  the  integral  (C.l)  have 
been  made  by  Hasselmann  (1963b)  and  more  recently 
by  Sell  and  Hasselmann  (1972).  The  latter  computa- 
tions, typically  30  minutes  per  spectrum  on  a 
CDC6600  computer,  lead  to  interpretations  of  the 
JONSWAP  measurements  that  strongly  suggest  that 
the  nonlinear  transfers  of  energy  play  an  essential 
role  in  the  development  of  the  wave  spectrum,  par- 
ticularly in  the  growth  of  the  wave  energy  at  low 
frequencies. 

Because  of  the  complexity  of  the  direct  calcula- 
tions of  (C.l),  Barnett  (1966),  using  a  Neumann  spec- 
trum, obtained  a  parameterization  of  (C.l)  as  part  of 
a  wave  prediction  scheme  based  on  the  earlier  Hassel- 
mann (1963a, b)  results.  Mitsuyasu  (1968)  applied 
this  parameterization  to  his  studies  on  nonlinear 
energy  transfer  in  wave  spectrum  and  obtained  fair 
agreement  between  calculated  and  observed  non- 
linear source  functions.  Resio  and  Vincent  (1976) 
satisfactorily  applied  Barnett's  prediction  scheme  for 
hindcasting  waves  in  the  Great  Lakes. 

Barnett's  parameterization  was  based  on  Hassel- 
mann's  (1962)  theoretical  results,  which  demon- 
strated that  the  wave-wave  interaction  of  a  spectral 
component  consists  of  an  active  part  that  transfers  its 
energy  to  other  components  and  a  passive  part  that 
receives  energy  from  other  components.  Thus, 
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Figure  C.l. — A  white  spectrum  and  its  computed  source  functions. 


with 


n        4.4  X  108EV       4,_     ../,     rt.,/o\3 
G=    \ —  cos   (0-  0O)(1  -0.42^  1 


(C.3) 


for/>  0.42  /0,  and  |0-0O|  < 

G  =  0,  otherwise; 

and 
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&M  =G-DS(f,9) 


(C.2) 


for/>  0.53  f0; 
D  =  0,  otherwise, 


1  +  16|cos(0-0o)|](/-O.53/o)3, 
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Figure    C.2. — Pierson-Moskowitz    spectrum    and    its    computed 
source  functions. 
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Figure    C.3. — JONSWAP    spectrum    and    its    computed    source 
functions. 


where  E,  f0,  and  90  are  total  energy,  mean  frequency, 
and  mean  direction,  respectively,  defined  by 


and 


E  =  \\S(f,d)dfdd, 

f0=h\fS{f.e)dfde, 


es(f,d)dfdd. 


(C.5) 
(C.6) 

(C.7) 


In  the  above  equations,  0  is  the  wave  direction  and  f 
=  w/2ir  is  the  linear  wave  frequency. 

We  applied  the  above  parameterization  by 
further  assuming  a  directional  spreading  factor  of 
(8/37r)  cos40  to  compare  some  of  the  results  given  in 
Sell  and  Hasselmann  (1972)  that  led  to  the  major  con- 
clusions based  on  JONSWAP  measurements.  The  re- 
sults for  a  fairly  white  spectrum,  a  Pierson-Mosko- 
witz spectrum,  and  a  mean  JONSWAP  spectrum  are 


presented  in  figures  C.l,  C.2,  and  C.3,  respectively. 
The  JONSWAP  study  (Hasselmann  et  al.,  1973)  con- 
cluded that,  as  the  sharpness  increases  from  figure 
C.l  to  figure  C.3,  the  evolution  of  a  sharp  peak  is  pri- 
marily controlled  by  the  nonlinear  energy  transfer 
evidenced  by  the  shifting  of  the  positive  lobe  in  Sell 
and  Hasselmann's  calculations  toward  lower  fre- 
quencies. While  not  necessarily  profound,  the  same 
conclusion  can  also  be  drawn  from  the  shifting  of  the 
positive  lobe  as  a  result  of  Barnett's  parameteriza- 
tion. Therefore,  with  an  overwhelming  savings  in 
effort  and  computer  time,  Barnett's  parameterization 
can  be  used  approximately  for  the  examination  of 
theoretical  nonlinear  energy  transfers,  particularly 
with  respect  to  the  location  of  the  positive  lobes  of 
the  parameterized  nonlinear  source  function. 
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Appendix  D:  Trispectra  During  the  Episodes  30  September  1972  and  7  October  1972 

We  present  here  the  calculated  trispectra  for  the         tral  sheets  located  30  minutes  apart,  corresponding  in 
episodes  of  30  September  and  7  October  as  figures         time  to  the  bispectra  shown  in  figures  16  and  17. 
D.l  and  D.2.  Each  figure  contains  six  sets  of  trispec- 
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Figure  D.la. — Trispectra  during  the  episode  of  30  September  1972,   2.5  minutes  from  beginning  of 
episode. 
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Figure  D.lb. — Trispectra  during  the  episode  of  30  September  1972,  32.5  minutes  from  beginning  of 
episode. 
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Figure  D.lc. — Trispectra  during  the  episode  of  30  September  1972,  62.5  minutes  from  beginning  of 
episode. 
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Figure  D.ld. — Trispectra  during  the  episode  of  30  September  1972,  92.5  minutes  from  beginning  of 
episode. 
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Figure  D.le. — Trispectra  during  the  episode  of  30  September  1972,  122.5  minutes  from  beginning  of 
episode. 
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Figure  D.lf. — Trispectra  during  the  episode  of  30  September  1972,  152.5  minutes  from  beginning  of 
episode. 
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Figure  D. 2a.  — Trispectra  during  the  episode  of  7  October  1972,  2.5  minutes  from  beginning  of  episode. 
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Figure  D.2b. — Trispectra  during  the  episode  of  7  October  1972,  32.5  minutes  from  beginning  of  episode. 
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Figure  D.2c— Trispectra  during  the  episode  of  7  October  1972,  62.5  minutes  from  beginning  of  episode. 
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Figure  D.2d. — Trispectra  during  the  episode  of  7  October  1972,  92.5  minutes  from  beginning  of  episode. 
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Figure  D.2e. — Trispectra  during  the  episode  of  7  October  1972,  122.5  minutes  from  beginning  of  episode. 
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Figure  D.2f. — Trispectra  during  the  episode  of  7  October  1972,  152.5  minutes  from  beginning  of  episode. 
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L  ABO 


The  mission  of  the  Environmental  Research  Laboratories  (ERL)  is  to  conduct  an  integrated  program  of  fundamental 
research,  related  technology  development,  and  services  to  improve  understanding  and  prediction  of  the  geo- 
physical environment  comprising  the  oceans  and  inland  waters,  the  lower  and  upper  atmosphere,  the  space  envi- 
ronment, and  the  Earth.  The  following  participate  in  the  ERL  missions: 


W/M         Weather  Modification  Program  Office.  Plans  and 
coordinates  ERL  weather  modification  projects 
for  precipitation  enhancement  and  severe  storms 
mitigation. 

NHEML  National  Hurricane  and  Experimental  Meteor- 
ology Laboratory.  Conducts  research  on  tropical 
meteorology  such  as  hurricanes  and  cumulus 
systems,  to  improve  understanding  and  fore- 
casting accuracy  and  to  develop  methods  for 
beneficial  modification  of  these  tropical  weather 
systems. 

RFC         Research  Facilities  Center.  Operates  instru- 
mented aircraft  for  environmental  research  pro- 
grams, and  provides  scientific  measurement 
tools,  logged  data,  and  associated  information  for 
airborne  and  selected  surface  research 
programs. 


WPL         Wave  Propagation  Laboratory.  Develops,  and 
applies  to  research  and  services,  new  methods 
for  remote  sensing  of  the  geophysical  envi- 
ronment. 


AL  Aeronomy  Laboratory.  Studies  the  physics, 

dynamics,  and  chemistry  of  the  stratosphere  and 
the  surrounding  upper  and  lower  atmosphere. 

NSSL      National  Severe  Storms  Laboratory.  Studies 
severe-storm  circulation  and  dynamics,  and 
develops  techniques  to  detect  and  predict  tor- 
nadoes, thunderstorms,  and  squall  lines. 


AOML     Atlantic  Oceanographic  and  Meteorological 

Laboratories.  Studies  the  physical,  chemical,  and 
geological  characteristics  and  processes  of  the 
ocean  waters,  the  sea  floor,  and  the  atmosphere 
above  the  ocean. 

PMEL      Pacific  Marine  Environmental  Laboratory.  Moni- 
tors and  predicts  the  physical  and  biochemical 
effects  of  natural  events  and  human  activities  on 
the  deep-ocean  and  coastal  marine  environments 
of  the  Pacific  region. 

APCL      Atmospheric  Physics  and  Chemistry  Laboratory. 
Conducts  research  on  mesoscale  meteorology 
related  to  atmospheric  pollution  and  inadvertent 
weather  modification,  and  to  the  prediction  of 
mesoscale  weather  phenomena.  Research 
subjects  include  cloud  and  precipitation  physics, 
diagnosis  and  numerical  models  of  mesoscale 
weather  events,  secular  changes  in  atmospheric 
chemical  and  particulate  composition,  and  radia- 
tive heat  transfer  applications. 

GLERL    Great  Lakes  Environmental  Research  Laboratory. 
Studies  hydrology,  waves,  currents,  lake  levels, 
biological  and  chemical  processes,  and  lake-air 
interaction  in  the  Great  Lakes  and  their  water- 
sheds; forecasts  lake  ice  conditions. 

ARL         Air  Resources  Laboratories.  Studies  the  diffusion, 
transport,  and  dissipation  of  atmospheric  pollu- 
tants; develops  methods  of  predicting  and  con- 
trolling atmospheric  pollution;  monitors  the  global 
physical  environment  to  detect  climatic  change. 


GFDL      Geophysical  Fluid  Dynamics  Laboratory.  Studies 
the  dynamics  of  geophysical  fluid  systems  (the 
atmosphere,  the  hydrosphere,  and  the  cryo- 
sphere)  through  theoretical  analysis  and  numer- 
ical simulation  using  powerful,  high-speed  digital 
computers. 


SEL         Space  Environment  Laboratory.  Studies  solar- 
terrestrial  physics  (interplanetary,  magneto- 
spheric,  and  ionospheric);  develops  techniques 
for  forecasting  solar  disturbances;  provides  real 
time  monitoring  and  forecasting  of  the  space 
environment. 
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